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Abstract 

We formulate and derive a generalization of an orthogonal rational-function basis for spectral expansions 
over the infinite or semi-infinite interval. The original functions, first presented by Wiener [30], are a 
mapping and weighting of the Fourier basis to the infinite interval. By identifying the Fourier series as a 
biorthogonal composition of Jacobi polynomials/functions, we are able to define generalized Fourier series' 
which, appropriately mapped to the whole real line and weighted, form a generalization of Wiener's basis 
functions. It is known that the original Wiener rational functions inherit sparse Galerkin matrices for 
differentiation, and can utilize the fast Fourier transform (FFT) for computation of the modal coefficients. 
We show that the generalized basis sets also have a sparse differentiation matrix and we discuss connection 
problems, which are necessary theoretical developments for application of the FFT. 

1 Introduction 

The approximation of a function by a finite sum of basis functions has long been a hallmark tool in numerical 
analysis. Over the finite interval much is known about expansion properties and periodic Fourier expansions 
or polynomial expansions are well-studied. On infinite intervals there are complications due to the unbounded 
domain on which approximation is necessary. Nevertheless many basis sets have been successfully investigated 
in this case; the Hermite functions provide a suitable method for approximation when it can be assumed 
that the function decays exponentially; for functions that do not decay exponentially, the so-called mapped 
Chebyshev rational functions can fill the void and open up the possibility for utilizing the fast Fourier 
transform (FFT); additionally, a Fourier basis mapped to the real line has been explored and provides an 
additional method for function approximation over the infinite interval. This last basis set serves as an 
inspiration for the family of basis sets proposed in this paper. 

Despite the available methods for function approximation over the infinite interval, there are shortcom- 
ings. The Hermite functions/polynomials do not admit an FFT exploitation and have problems approxi- 
mating functions that do not decay exponentially (which is to say, most functions). However, the solutions 
to differential equations have been relatively successful by Hermite approximation and in some cases give 
superior approximations when compared to a Chebyshev (mapped or truncated) approximation [5]. The 
Whittaker cardinal interpolant functions [29], or Sine functions, provide a remarkably simple method to 
approximate a function with known equispaced evaluations. The drawback is a relatively small class of 
functions for which such an expansion is complete. However, the ease of applying Sine methods has led to a 
great number of applications [5T] . The Chebyshev rational functions [I] , [TU] are robust with respect to the 
deficiencies of the Hermite and Sine bases, but they have some disadvantages compared with the generalized 
Wiener basis we will derive. 
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On the semi-infinite interval Laguerre polynomial/function expansions are the classical approximation 
technique [T], but these techniques suffer from the same problems as Hermite expansions. An alternative 
technique involves mapping Jacobi polynomials to the infinite interval [7J. This mapping technique makes it 
possible to accurately approximate algebraically-decaying functions on the semi-infinite interval, but intro- 
duces some computational issues for the solution to differential equations. The generalized Wiener basis can 
be employed on the semi-infinite interval; this results in a basis set that is also a mapped Jacobi polynomial 
methods. However, the Wiener mapping is very different from that presented in the literature, and therefore 
acts as a competitor to these existing techniques. 

Our generalized basis is inspired by a collection of orthogonal and complete functions originally proposed 
by Wiener [30] . He introduces the functions 

(1 - ix) n 

as Fourier transforms of the Laguerre functions. He furthermore shows that these functions are orthogonal 
under the L? conjugate inner product. Higgins |19] expands this result by presenting the functions <j) n 
along with their complex conjugates as a complete system in I? . Following this, others have followed up 
on these functions by applying them to the solution of differential equations [13], [11]. We note that the 
functions 4> n (%) presented above have magnitude that decays like ^ as |a;| — > oo. We will generalize the 
above functions so that they have decay ^ for any s > \. The ability to choose the rate of decay of the 
basis set is an advantage if such information is present about the nature of the function to be approximated 
or the differential equation to be solved (e.g. [20], [22]). Furthermore, we will show that this basis admits 
sparse Galerkin matrices and that the fast Fourier Transform can be used in certain cases to evaluate and 
manipulate the series. 

This paper is concerned with the derivation and theoretical properties of the generalized Wiener rational 
function basis. Computational considerations, numerical examples, and comparisons with existing basis sets 
are presented in a second part. The outline of this paper is as follows. In Section [2] we formulate and derive 
the basis, which is heavily based upon a generalized Fourier series. Section [3] follows with some properties of 
the basis functions based on their close relationship to the canonical Fourier basis, and Section [4] concerns 
the properties that can be derived from the relation to Jacobi polynomials. In Section [5] we discuss how 
the Wiener basis set may be used to approximate functions on the semi-infinite interval. Finally, we briefly 
present mapped Jacobi polynomials as an alternative method in Section [6] and summarize and present an 
outlook in Section [JJ for Part II, dealing with numerical issues. 



2 Derivation of the basis 



We begin by stating the major goals and the path we will take in accomplishing those goals. We seek a 
collection of L 2 -orthogonal and complete basis functions whose domain is the entire real fine. In addition, 
we desire the ability to specify a parameter s > \ that will denote the polynomial decay at ±oo of each of 
the the basis functions. 

Drawing inspiration from Wiener and his orthogonal basis functions, we seek a collection of functions 
d>[ s \x) for x £ R and k € 7L such that \d>\f \ is a complete, orthogonal system for any valid s. Our 

method relies on the observation that the functions |l]) are weighted maps of the canonical Fourier basis e mB 
for G [0, 2ir] (see e.g. [TU], [25]). We will first generalize the Fourier basis on [0, 2ir] so that it will have the 
properties we desire on the infinite interval; we will then map the generalized Fourier basis to the real line 
and weight it accordingly to achieve the desired rate of decay. 



2.1 Notation and setup 



We shall reserve the variables x, z, 9, and r as independent variables on certain domains and list the domains 
and transformations in Table [l] The variable r G [ — 1 , 1] is the standard interval over which the Jacobi 
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r = cos 9 
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x ~ v l+r 


^ gi arccos r 


= arccos r 


re [-1,1] 



Table 1: Isomorphic transforms between different domains. 

polynomials are defined. The interval 9 € [0,7r] is the image of the r interval under the map 9 = arccos r. 
The variable z S T + is the upper-half of the unit circle in the complex plane, and x 6 [0, oo] is the positive 
half of the extended real line. 

In much of what follows we will mix notation and write expressions both in terms of e.g. r and 9. It 
should then be understood that r = r(9) and/or 9 = 9(r). Furthermore, we will extend the domains of 9, z, 
and x to be [— tt, tt], T, and R, respectively later in the paper. 

We denote L 2 (A, B; w) = L 2 W (A, B) the space of square integrable functions / : A — > B under the weight 
to. We endow L 2 ^ (A, B) with the conjugate bilinear inner product; the notation for this inner product is 
(•, ■,} w - The omission of w indicates the unit weight measure. The norm on this space will be denoted |j -|| . 
The following weight functions will be used extensively in this article: 



w ( a >P)(r) = (l- r ) a (l + rf 

w£' 5) (d) =t»(W(r(9)) = {I + cos 9y {I- cos 9) s 



w^\x)=w[ s ' t \9(x)) = - 



■ta+t 



■21 



(1 +x 2 ) s V(l +x 2 ) 

In addition, we will make use of a phase-shifted square root of Wx and wi' t ' S \ which we define as: 



(s,i) 

w x exp 



i(s + t) 



(n-9(x)) 



{]w^ S) {9) = yw£' S) (x(0)) = 2(^) si 
2.2 Jacobi polynomials 



sm d { "- \ cos 1 I — | exp 



(x - i) s+t 



(tt-0) 



(2) 
(3) 



The classical Jacobi polynomials P« a '^ are a family of orthogonal polynomials [27] that have been used 
extensively in many applications due to their ability to approximate general classes of functions. They are 
a class of polynomials that encompass the Chebyshev, Legendre, and Gegenbauer/ultraspheric polynomials. 
These polynomials will form the building blocks for our generalization. 



The Jacobi differential equation is 

(1 - r 2 )p" + [0-a-(a + /3+ 2)r] ft + n(n + a + [3 + l)p = 0, re [-1, 1], 



(4) 
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and for a,/3 > — 1, n £ N the only polynomial solution p = Pn°'^\x) is a polynomial of degree n. 
The restriction a,f3 > —1 is necessary to ensure integrability of the weight and thus existence of an 

L 2 -constant function solution. The family of polynomials { P^^\x) \ is complete and orthogonal in 



L 2 (\-\,\],Wi\wi a ^. We denote hfi 



(«,/3) 



P 



(«,/3) 



and define the normalized polynomials as 



The orthonormal Jacobi polynomials P^ 01 '^ will be integral in the derivation of the Wiener rational function 
basis on the real line. In addition, we require a minor generalization of Jacobi polynomials: we perform a 
change of the dependent variable in Q to obtain: 

Lemma 1 (Jacobi Functions) The Jacobi functions defined as 

p(a,p,a,b)^ = (1 _ r) a (1 + ^bp(a,P)^ 

satisfy the following properties: 

1- \ Pn a '^' a ' b \r) \ are orthogonal and complete in L 2 ([—1,1], R; 2a ^ 2fc M. 

L J nSN V / 

2. The Prf" , P' a "' b \r) are eigenfunctions p n (r) of the Sturm- Liouville problem 

d 



\p(r)p'(r)] + q{r)p(r) - \ n w(r)p(r) = 0, 

Ul 

which is defined by the parameters 

p{r) = (1 - r ) a+1 - 2a (l + r f +1 - 2h 

q(r) = [a{a-a){\-r)- 2 + b{(3-b){\ + r)- 2 ]{\~r) a+l - 2a {\ + rY+ 1 - 2b 
w(r) = (l-r) a - 2a (l + r)P- 2b 

X n = n(n + a + (3+ 1) - 2ab+ a(/3 + 1) + b(a + 1) 



(5) 



(G) 



The proof is mathematically simple but algebraically tedious and we omit it. We shall actually only require 
the result of Lemma 1 for a = b = g. Many of the results in this paper require the use of numerous recurrence 



relations involving Jacobi polynomials; these relations are given in Appendix [A| equations (26l-(33l. 

The idea behind the formation of the Jacobi functions introduced in Lemma [l] is not novel and has 
already found use in the literature. In |17] the 'generalized Jacobi polynomials/functions' are denoted jn' , 
and are defined for all a, (3 £ K as 



p (-a -0,-a,-0) 
p(-a,/3,-a,0) 

Til > 

p ( a -f3,O,-0) 
p(a,/3,0,0) 

TL 1 , 



a < -1 and (3 < -1, 
a < — 1 and (3 > -1, 
a > -1 and (3 < -1, 
else, 



where the index ni is defined as 
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n-[-a\ - L-#], 
n - [-a], 
n- L-/3J, 



a < -1 and /? < -1, 
a < -1 and /? > -1, 
a > — 1 and /? < — 1, 
else, 



and the integer floor function is denoted [-J . These functions are only defined for certain values of n but 
|17j presents significant approximation theory using them. They are advantageous for solving high-order 
differential equations with boundary conditions via a global spectral expansion. 

Finally, we present two classical notational conventions that we will use briefly in the next section. The 
classical Jacobi polynomials that result from the cases a = {3 = — | and a = /3 = +\ are the Chebyshev 
polynomials of the first and second kinds, respectively. Recalling the relation r — cos 9, these polynomials are 
typically denoted T n (r) and U n (r) and they have a very special and concise representation as trigonometric 
polynomials: 

%P^T 1 ^ 2, ^ 1 ^' 2 \r) = T n (r) = cos(n9) = cos [n arccos(r)] 



7T p(l/2,l/2) 
2 r " 



(r) 



U n (r) 



sin[(n+l)0] 
sin 



sin[(n+l) arccos(r)] 
sin[arccos(r)] 



2.3 Generalizing the Fourier basis 



In this section we will generalize the canonical Fourier basis given by 

= e tke . 

Our methodology is based upon the following dissection of the Fourier basis for k ^ 0: 

cos (k9) + ism(k9) 



ik-0 



cos(|fc|0) + isgn(fc)siu(|fc|0) 

Ij fe |(cos0) + isgn(fe)sin(0)Z7| fe |_ 1 (cos0) 



p(-l/2,-l/2) 
^\k\ 



(COS0) 



* sgn(yfc) sin(0)^ ( ? ! | / J ! ; 1/2) (cos ( 



(a) 



('<) 



We have broken down the Fourier basis into two components: the first component (a) is even with respect 
to 9 as it is simply a polynomial in cos0. The second term (b) is odd in 9 as it is a polynomial in cos0 (an 
even function) multiplied by the odd function sin0. This breakdown suggests that we can construct more 
general kinds of Fourier-type functions by augmenting the type of polynomials employed. 

However, we cannot switch around polynomials with impunity; we want to retain orthogonality (at 
least with respect to some weight function). The separation into terms (a) and (6) above elucidates the 
biorthogonal decomposition of the Fourier basis. The (a) functions are orthogonal with respect to each 
other, and with respect to the (6) functions. In this case, the biorthogonality is manifested as an even-odd 
separation. Suppose we wish to generate a basis set orthogonal under the weight 1 + cos 9 = l + r. Naturally 
we can do this for basis (a) by changing the second Jacobi class parameter from (3 = — | to /3 = +|. In 
order to do this for basis (6), we use Lemma [T] 

For a,@ > —1, we have the polynomials P^C'^ that are orthogonal in L 2 ([— 1, 1],R; w^^'j . By setting 



we also observe that the Jacobi functions p ? ( Q+1 '' 3+1 ' 1 / 2 ' 1 / 2 ) 



a = b = ^ in Lemma 

are orthogonal under the same weight. If we set a 



(l-r 2 ) 1 ^' 



(a+1,/3+1) 



4, and add these two functions together with the 
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appropriate scaling factors, then we exactly recover the Fourier basis by reversing the dissection steps above 
(i.e. by creating a biorthogonal construction). Of course, we are free to choose any values of (a, (3) that we 
desire in order to derive generalized trigonometric Fourier functions. In fact, this technique has already been 
used by Szego f27 to determine orthogonal polynomials on the unit disk. Because the statement in [27] is 
merely a passing comment and is a markedly different result than what we desire, we present the following 
theorem: 

Theorem 1 (cf. Szego, [27 ) For any 7 > — |, the functions 

i o (-iA7-Va) (co8(?)i fc = 

(7) 




P 



1 (cos 9) + i sgn(fc) sin(0)P | ( fc 1 | / _ 2 ' 7+1/2) (cos 9) 



fc^O 



are complete and orthonormal in L 2 ([— tt, tt], (D; Wg 1 ' " 1 
Proof For orthonormality, it suffices to show 

1 / n(-l/2, 7-1/2) , Q s 5(-l/2,7-l/2) , a \\ 

1- (P^ J (cos0),P^ ; (cosfl)^ =2«| fc |,|i| 

2. (sineP\l(^ +1/2) (cos9),sin9P^ +1/2) (cos0))^ Q) = 2<5| fcMi |, for k,l + 0. 

3. ^AlC-V^T-i/a) (cosfl) .sineP^'^^fcoBe)) _ 0) - 0, for I ± 0. 

The first property is a direct result of orthonormality of the normalized Jacobi polynomials P and the 
observation that on [0,n], (/(cos0),<7(cos0)) ( 7 ,o) — (/(?"), <?(r)) (-1/2,^-1/2). The second property is a result 
of the same observations as the first property along with the result of Lemma [l] The third property results 
from the fact that an odd function integrated over a symmetric interval is 0. Orthonormality then follows 

from an explicit calculation of /^/ < ^ / \^/ < f " > \ using the above three properties. 

For completeness we note that any function / 6 L 2 can be decomposed into an even f e and an odd f a 
part. That f e is representable is clear from the fact that P„ 1 ^ 2n 1 / 2 - 1 (cos 6) is complete over 6 6 [0,7r], 
which by symmetry implies completeness over all i 2 -even functions fe. Similary, the collection of functions 
sin 6Pi~ 1/2 ' 7 ~ 1/2) is complete over all X 2 -odd functions f a by Lemma 111 Linearity and orthogonality of the 
even and odd parts yields the result. □ 



Remark 1 Szego [27 gives a more general result that involves orthogonality over the weig ht W y' d) for 5^0. 
We do not require this level of generality; for 6^0 the weight function becomes zero at 9 = 0, which we will 
see does not help our cause. Indeed, it is possible to generalize Szego's result: he derived polynomials on 
the unit disk orthogonal with respect to w e 7 ' S \ By using Lemma^with a, b different from |, we can derive 
non-polynomial basis sets that are orthogonal under a great variety of weights. These functions naturally 
may not be periodic on 9 £ [— tt, tt] if the quantity (1 — r) a (l + r) b cannot be periodically extended in 0-space 

to [— 7T, n]. 

We will refer to the functions as either the generalized Fourier series, or the Szego- Fourier functions. In 
the definition of the functions Vt^ it is desirable to use the L 2 -normalized versions of the Jacobi polynomials 
P, rather than the standard polynomials P. If the standard polynomials are used, then the norm of the 
Szego- Fourier functions depends on the rather unpleasant-looking sum h^^ 2 ' 1 1 ^ 2 ^ + h^ 2 '^ +1 ^ 2 \ and 

using this convention implies that 'I'j^ is not orthogonal to ^^l- 

We can also distribute the weight function onto the basis functions, which gives us orthogonality in the 
unweighted L 2 -norm: 
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Figure 1: Plots of the weighted Szego-Fourier functions ipj?'(0) for k = 0,1, 2, 3, and 4. Real part (top) and imaginary part 
(bottom). 



Corollary 1 For any 7 > — ^, the functions 



V w e"" p(-l/2,7-l/2) 
x/2 *0 



(cos 0), 



(-1/2,7-1/2) 
|fc| 



(cos 9) + i sgn(fc) sin(6»)P | l fc y_ • 



(1/2,7+1/2) 



(cosfl) 



k = 



fc^O 



are complete and orthonormal in L 2 ([— 77, 7r], <D). 



Due to the properties of y Wg 7 ' '' given in (3 1, the functions ij}^\ff) decay like (cos |) 7 at 9 = ±tt. This is 
exemplified in Figure [l] where we plot the real and imaginary parts of the functions for 7 = 2. The even/odd 
behavior in 9 for real/imaginary components depicted in the figure depends on the even/odd parity of 7. 
(There is no such characterization possible when 7 g" N .) Clearly for 7 = we have vD^ — tp^ — y|^e lfc6 ', 
the canonical Fourier basis. 



2.4 Mapping to the real line 

Having developed the necessary preliminaries on the finite interval, we now jump to the infintc line x £ HI 
using the mappings introduced in Table [T| To facilitate the mapping, the following identities characterizing 
the mapping between (9-space and a;-space are useful: 



cos 9 



sin 9 



l-x 2 
l + x 2 



2x 



(1 — cos? 
(1 + cosf 



2x 2 



7 



Using these identities, we rewrite and relabel the functions ^^(9): 



*i'- X) (fi(x)) 



j_ p (- 1/2,8-3/2) ( i- x < 



P, 



(-1/2,8-3/2) / l_x 



|fc| 



l+z 5 



2ixsgn(fc) p( 



,(1/2,8-1/2) ( X -x< 

|*|-i v i+3:2 



fc = 



The above definition is valid for any s > |. s = 1 corresponds to a mapping of the canonical Fourier 

basis (i.e., s == 7 + 1). These functions are orthogonal over the weight w x s '°\ By following the route from 
Corollary [T] we can distribute the weight over the basis functions, and in this particular instance we choose 
the phase-shifted square root given in 



w 



2 (^\ £(-1/2,8-3/2) (W\ /, = () 



(l+"0 

f l-x 2 \ , 2«sgn(fe) 5(1/2,8-1/2) f l-x 2 \1 h , n 

^1+k 2 ^ + x 2 + l ^|fe|-l \^l+z 2 ,IJ' ' U - 



{x-iY 1 

(8) 

2(i~ 1 ) [p(-l/2, 8-3/2) 

The functions (JsJ) are what we call the generalized Wiener rational functions. At present there is no clear 

reason why we have chosen to use \J w x s '°^ instead of the usual square root \J w x s ^ to distribute the weight. 
However, the corollary following the coming proposition should provide part of the motivation. 

Proposition 1 For any s > |, the functions <$>*£\x) are complete and orthonormal in 1? ^R, (D; wi 3 ' ^ 

The functions (x) are complete and orthonormal in L 2 (R, (D). Furthermore, the decay rate of these 
functions can be characterized as 



lim 

\x\ — »oc 



x l ^\x) 



< 00, t < s 



Corollary 2 Recalling the definition of Wiener's original basis functions (/> n (x) in the following relation 
holds: 

iV2(f>^(x) = 4> n (x), n€N . 

We show plots of the functions (f>^ in Figure [^J The conclusion of the corollary is easily seen if one makes 
the connection 

%9 _ 1 ~ x 
1 + x 

along with knowledge of the fact that 3?^ (#) = ipf\9) = ~^^G lkB . We have thus shown that the orthogonal 

functions <j>^ over the real line are a generalization of Wiener's original basis set. Furthermore, <f>£ decays 
like x~ s while retaining orthogonality under the same unit weight measure. When s is an integer, the 
functions are also purely rational: they are the division of one complex-valued polynomial in x by another. 
This connection was rather helpful in the nascent stages of the computing when the calculation of a non- 
polynomial function required significantly more computational investment, but now this property is probably 

more aesthetic than functional. As a result, our use of the quantity y 1/^'°^ is not entirely necessary for 



purposes of evaluating the functions; it is equally valid to use the traditional squre root \l w x s '°\ 

By using the traditional square root, one sacrifice made is that the analogous written form of Corollary [2] 
becomes less fortuitous and is complicated by .T-dcpcndcnt phase-shift factors. The same observation is true 
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X 



Figure 2: Plots of the functions 0^ 4) (x) for k = 0, 1, 2, 3, 4. 



of the weight y Wg 7 ' ^ used in the definition of the Szego- Fourier functions ^^{0) in Corollary 



A second 



reason to use the phase-shifted square root is that it can be written in the following convenient Torm: 

(9) 



0,o) 



The utility of this expression will become clear when we consider the connection problems in Section [4] 

We have accomplished our goal of deriving basis functions satisfying tunable decay rate while maintaining 
L 2 -orthogonality. However, it is not clear that these are superior or useful functions. We will now present 
some properties of the basis and make the argument that these basis functions indeed are very useful for 
solving problems in scientific computing. 

3 Fourier-Derived Basis Properties 

In this section we explore some of the desirable properties of the generalized Wiener basis set ^<p^ j , 

s > | based on their close relation to Fourier Series. The argument we make is that these functions inherit all 
the useful properties of the Fourier basis with the additional property that the decay rate s at |ar| = oo may 
be chosen. Many of these properties (e.g. the sparse modal differentiation matrix) rely on Jacobi polynomial 
properties covered in the next section. In particular, although application of the FFT is indeed a virtue of 
this basis, we will discuss it only in Part II, which focuses with computational issues. 

3.1 Symmetry 

The derivation of the basis functions automatically yields various simple properties. Note that due to the 
mapping, any property of the basis on the real line x £ R also applies to the respective trigonometric interval 
9 6 [— 7r, 7r]. We omit the proof of these properties as they are elementary: 
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1. Index symmetry 



= (10) 



1^)1 - 1^(^)1 



<^(x) = £_jU(*) 

2. Function symmetry 



Re 
Im 



{^\x)} = He{*W(-x)} 

{4 S) W} = -Im{<f>i s) (-x)} 

I^ S) WI = l*i 8) (-^)l 

I^ S) WI = I^ S) (-^)I 



3.2 Periodicity 



Because trigonometric polynomials are periodic over 8 £ [— 7T, 7r], we cannot expect this condition to be 
violated on the infinite interval x £ R. From the viewpoint of expanding functions over the infinite interval 
R, the points x = ±co are both unique points. However, because of the mapping, the basis functions view 
the points x — ±oo the same as they view the points 8 = ±7r: i.e. they are the same point. This serves 
as a disadvantage if we wish to e.g. expand functions with different decay rates at ±oo because this is in 
effect non-smooth behavior of the function at a single point, which degrades the convergence rate of the 
approximation . 

In particular it is known that although a Fourier series approximation will converge in the L 2 sense for an 
L 2 function, the rate of convergence is only algebraic if the function is non-periodic. Naturally, this deficiency 
will follow us to the infinite interval. Indeed, such observations have already been made |12j . Empirical 
studies we have carried out show that the concern of periodicity is not paramount and frequently one can 
overlook it when comparing results to other basis expansions. Nonperiodic behavior is often manifested as 
algebraic decay at x = ±oo, where existing basis sets already have problems in approximation. In Part II 
we will present examples that explicitly illustrate this lack of fast convergence rate when the function to be 
expanded is not 'periodic' at x = ±oo. 

Note that although it may seem a bit unnatural that periodicity is a condition at x = ±oo, in fact it is 
not surprising at all. One mayconsider our mapping as a rather unremarkable tangent map from 0-space 
to x-space as written in Table 111 However, it is more deep than that: the functions ^^(6) and tp^id) are 
periodic basis sets for complex- valued functions on the unit circle T. In other words, we may actually view 
these basis sets as functions of z £ (D . What looks like a tangent mapping from #-space to x-space is actually 
a linear fractional map (a Mobius transformation) from the unit circle (the complex plane) to complexified 
x-space (the complex plane). 

Linear fractional maps are structure-preserving maps of the complex plane; an illuminating way to 
consider this is by identifying the complex plane (D with the Riemann Sphere (see Figure |3(a) I . Then the 



linear fractional map we've chosen merely takes one great circle (the unit circle in z-space) to another great 
circle (the real line in complexified x-space). Therefore, our approximation is nothing more than a rotation 
of functions on the Riemann Sphere (see Figure |3(b)[ ) ; the target space simply happens to correspond to the 
real line. From this point of view, periodicity at \x\ — oo (analyticity at complexified x = oo) is natural. 
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(a) Illustration of the stereographic connection between the Riemann Sphere and the 
complex plane. The equator corresponds to the unit circle, and the meridian can be 
identified with the real (x) axis. 




(b) The effect of the linear fraction map we've chosen to take 8 = arg^i to x = Re{z2j: a rotation 
of the Riemann Sphere. 

Figure 3: The linear fractional mapping that relates x to 8 has an illuminating representation when viewed as a transformation 
of the complex plane. 



Nevertheless, this 'natural' periodicity can be problematic if we attempt to approximate a function that 
is not complex- analytic at x = oo. In Part II we present examples of functions that are not analytic at 
x = oo and we will empirically analyze the impact of violating the assumption of periodicity. 

4 Jacobi-Derived Basis Properties 

The generalized Wiener functions are composed of Jacobi polynomials, and so it is reasonable to expect that 
we can use the properties of the Jacobi polynomials to perform certain tasks using the Wiener basis. Indeed, 
we can form recurrence relations, connection coefficients, a Gauss-like quadrature, and obtain an extremely 
useful sparsity result for the Galerkin stiffness matrix. 
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4.1 Recurrence Relations 



Due to the strong dependence of the Szego- Fourier functions on the Jacobi polynomials, they inherit six-term 
recurrence relations from the three-term recurrences for orthogonal polynomials. 



*n+l 



1 n 



B 



(7) 



(7) r /'-i (7) ,7,(7) j- /^Mtf/M 

°" n — 1 — rt — (n— 1)' 



(7) 
n+1 



[/„ 7) cos( 



(7) 



(7) 



'crW cos e-viZ 5 



(7) 



(7), T ,(7) 



K -n -(n-1)' 



(7) 
n+1 



t/„ 7) zsin0-K (7) 



(7) 



^(7) + 



r(7), T ,(7) 



K -n -(n-1)' 



We give formulae for all the real- valued constants A,B,C,D,U,V,W,U ,V ,W in Appendix |a| Note that 
since the *S>^ are not polynomials in z = e t6 , there is not a three-term recurrence as there would normally 
be for orthogonal polynomials on the unit disk (unless of course 7 = 0). Although the above formulae are 
complex- valued six-term recurrence relations, they are no more difficult computationally than the pair of 

three-term recurrences necessary to generate Pn"'^ and p^ a+1 '^ +1 ) because \E'„ 7 ' ) is the complex conjugate 
,(7) 



of V?_„ and therefore does not need to be generated independently. Direct use of any of the above six 

,(7) 

- fc 



term recurrences for generating the is just as expensive as forming tyY" by the even/odd synthesis of 



Pn '^ and p( a+1 <P + ^ m Theorem [l] However, using presumably existing routines for evaluating Jacobi 
polynomials and then synthesizing them is likely easier from an implementation view. 



It is reassuring to note that simplifying the recurrence constants in the case 7 = yields, up to normal- 

(0), 

k 



/2tt 



ization, the trivial recurrence relations for the monomials on the unit disk ^^(arg z) 



(o) 

n+1 



(o) 

n+1 



2co fl M#>-*W i) 



Naturally, a recurrence relation for the unweighted \&1 7 (6) translates directly into one for the unweighted 

Wiener rational functions $^\x). The weighted functions ip^\o) and <fi^\x) can be generated by first 
generating the unweighted functions and then multiplying by the phase-shifted square root f/w. 



4.2 Connection Problems 

One advantage in using the generalized Wiener rational function basis is the ability to choose the parameter 
s, which indicates the rate of decay. In many applications, it may be useful to augment the basis functions 
mid-computation to suit the dynamics occuring at a particular time. In this case, one would like to be 
able to transfer from one basis to another while keeping the (finite-term) function expansion identical. 
We will also see in Part II that this problem also appears in an algorithm utilizing the FFT. In classical 
orthogonal polynomial theory, the problem of equating one expansion to another boils down to determining 
the connection coefficients. Before undertaking this task, we first outline the major tasks we wish to perform. 

There are two main tasks on the infinite interval that require connections of some form: 

1. Usage of the fast Fourier transform - transforming N nodal evaluations into N modal coefficients (or 
vice-versa) for an expansion in cf>( s \ 

2. For a given expansion in c/)^ (i.e. a set of modal coefficients), translating this into a modal coefficient 
expansion in (j)^ for some s ^ S. 
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In Part II where we outline computational considerations, we will address the above tasks. However, for now 
it suffices to note that these two tasks can be reduced to the following three connection problems in 0-space: 



1. The vI/(T)-vlK r ) connection (a necessary ingredient for all connection-like tasks) 

2. The ^("-^(t) connection (a generalization of the FFT task) 

3. The i()^-tp^ r ' connection (identical to modification of s) 



In Sections 4.2.1j|4.23 we will tackle each of these problems. Note that modification of any of the following 



finite-interval algorithms for the infinite interval is trivial: the relations ^f^\0) = Q 1 ^ ( x )> ^ki®) = 
^(z)' an( i 7 := s — 1 allows for us to easily employ the same operations, whether we want to do it in 
6>-space or a;-space. 



4.2.1 The Connection Problem 



Suppose we have a function / £ L 2 f[— n, it], (D; w^^j f]L 2 h— n, n], <D; wf*^ with a Fourier expansion for 



some 7 > — I 



The goal is determine a way to re-expand / in a Fourier expansion for a different decay parameter T 

;( r ), T ,(n 

fee 



The shift T — 7 can take values in the interval (— | — 7, 00). Naturally one may equate the two expansions 
and use orthogonality to relate one set of expansion coefficients to the other: 



/f ) = Er ^ U (r) 
We can then define the connection coefficients 

where we have suppressed the dependence of A on 7 and T. Our task is to determine how to calculate these 
connection coefficients. Due to orthogonality, it is clear that 

A^=0, |J|<|fc|. (11) 

This implies that the connection problem is solved via the relation 

l e z, 
|/| > |fe| 



Relation ( 12 1 is still not attractive: we must perform an infinite number of operations for an exact connection. 
If we only have a finite expansion (say a total of TV modal coefficients) , we must still perform 0(N 2 ) operations 
to capture all the information at our disposal. However we will show that, for integer values of the shift 
r — 7, the connection problem can be solved inexpensively. To be precise, we will show that for G G N, ( 12 1 
reduces to 

fl 1+G) = E fi J)x tv (1 3 ) 

fe+G>|;|>|fc| 
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That is, only 2{G + 1) operations per coefficient are necessary to solve the connection problem (independent 
of k, and of any truncation size N). We refer to the above collapse of the infinite connection problem (12 1 



into the finite iV-indepedent problem (13 1 as a sparse connection 



In order to relate one Fourier function to another, we first recall a result from |25j using (29 1 - (32 1 that 



states that the connection coefficients binding one Jacobi polynomial class to another are sparse in certain 
special circumstances. 

Lemma 2 For any a, (3 > —1 and any A,B,€ N , the connection problem 

oo oo 

f(r) = £ ft^Pt 0) (r) — f(r) = £ f^ A ^P^+ B \r\ 



n=0 



n=0 



can be solved exactly via the relation 



f, 



(a+A,(3+B) 



= E A 



(14) 



In the above we have suppressed the dependence of A p on a, f3, A, and B, but in the sequel we shall occasion- 
ally refer to the above coefficients as A„ m (a, /3, A, B). The result (14) is not a trivial one; the upper limit 
for the sum on the right-hand side is oo for a general connection problem. For the very special cases satis- 
fying the lemma, the exact connection becomes finite. We have not shown how to obtain the Jacobi-Jacobi 
connection coefficients X p . For this, one may use explicit formulae given in |23j or [2], or one may utilize the 
algorithm given in 25\. 

The above result can be expanded to apply to the Szego- Fourier functions ^ (0) and the corresponding 
mapped functions $^\x). 

Proposition 2 For any 7 > — | and any G E IN, the connection problem 



2 

00 



m= E I 



(7+G)^,(7+G) 



(*), 



k=— 00 



k= — oo 



can be solved exactly via the relation 



\k\+G -\k\ 

;r G) = e a m// 7) + e a m// 7) - 

2=|*l l=-\k\-G 



(15) 



Note that |l5| is exactly (13 1. By making the connection s — 1 < — > 7, we recover A*j = A*j, where &^\x) 

are the maps of the Szego- Fourier functions . We stress again that this result is nontrivial. This also 
yields the functional connection 



< G 



y 



G<\k\<7 



^, m ^ +G) (0), \m\>G, 



(16) 



(7) , 



is a linear combination of at most 2G + 1 functions * 



(7+G) 



Note that the Fourier relation (16 1 



parallels (151 in exactly the same way that the Jacobi relations (29 1 - (30 1 parallel (31 1 - (32) 



We now illustrate how to calculate the Szego- Fourier connection coefficients A* in Proposition [2] from 
the Jacobi coefficients A p . In the following, we make use of the notation: 



1*1-1. 



P = l 
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Figure 4: Illustration of steps taken to perform ty-ty connections. The operator += is the addition-assignment operator. 



From the definition of \E'1 7 ' > in (|7l) we have 



(«,/3) _ ^(7) 



(a+1,/3+1) 



(7) 

-k 



n>0, 



Ifcl -Ifel 

P (a,/3) = 



Therefore, from the modes we can derive two sets of Jacobi modes: 



Ja,0) _ ?( 7 ) , ?(7) 
e n — Jn 1 J-n , 



n > 1, 



.(a+1,/3+1) 



i (7) _ ?(7) > Q 

in+l 7-n-l' " — u ' 



= V2ti 



(7) 



The Jacobi modes e„ are modes in an expansion in polynomials P„ ^ and the modes b n are for an expansion 
in pj 7 a+1 'f 3+1 \ with these mod 
coefficients using Proposition [2] 



in p^ Q+1,/3+1 \ With these modes in hand, we can use the Jacobi connection coefficients to promote the 



.(a,/3+G) 



,(a+l,/3+G+l) 



— Sm=0 ^ 



E™=o A^„ + ™o I ( l Q + ^' /3+lj , where A p = A p (a + 1, /3 + 1, 0, G), 



where A p = A p (a, /?, 0, G), 



n > 
n > 0. 



Finally we redistribute the modes back into Szego- Fourier form to yield what we desired: 



f(7+G) 
Jn 



f(7+G) 
J —n 



(a,/3+G) ~(a+l,/3+l+G) 



s (a,/3+G) ~(a+l,/3+l+G) 



n > 1 

71 > 1 



Jo 



(7+G) 



~( = ,/J+G) 



The whole procedure is illustrated graphically in Figure [4] We may explicitly write the connections as: 
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(7+G) 



G 
m=0 



2 [ Ifcl.lfc 



0) 

\+m 



(7) 

|fe|+m 



G 
m=0 



P;(a,/3) 
|fe|,|fe|+m 



s f rr n ('i-n P;(Q+1 ' /3+1) 
SglHfcJA |feUfc|+m 



?(7) 

-|fc|-m> 



Ifcl > 1 



G , G 



fh+G) _ J_ v-^ x p K",/3) f(7) , lr fh) 

JO ~ /o 0,m Jm ~< r~ Z_^i A 0,m J-m" 

m— m=0 

Therefore we have an explicit expression for the Szego- Fourier connection coefficients in ([2]): 



A P;(a,/3) ±sm m A ^(«+l,/3+l) 1 1. 1 > 1 

A lfcl.lfcl+m ±& g n W A |fc|.|fc| +m i l K l ^ 1 



l |fe|,|fe|+m - 1 - °5"W^|fc|,|fc|+ m 
A k.±flfcl4-m1 = \ ( 17 ) 



J,P;(a,/3) 



Of course, owing to observation ( 11 ), the above equation restricts < m < G. As mentioned, this connection 



relation is also valid for converting an expansion in the functions <$> k (x) to one in the functions $^ +S \x) for 
SsN since the modes for these two expansions are the same. Let f(0) be given and define g(x) = f(8(x)). 
Then for all 7 > § : 

This completes the connection problem. The reverse connection problem (converting fu + modes 

to modes) is solved by reversing the above procedure (all steps are invertible) and use of the fact 

that the forward Jacobi connection problem with integral separation is banded upper-triangular and thus 
the backward connection is O(N) calculable sequentially via back-substitution. See [25] . 

We have determined how to quickly and exactly accomplish the connection problems for the unweighted 
functions 

E fc ^' (s) *i s) (x) — E k st Xs+s) H s+s \^ 

in O(N) time where N is the total number of modes when S,G G %. These connections can be performed 
by utilizing the connection coefficients in ( 17 1 along with the sparse connection result of Proposition [2] For 
S,G *E, there is no sparse connection result for the modes, and so while the connection coefficients A* ; 
can still be calculated based on known connection coefficients for Jacobi polynomials, the coefficients do not 
terminate finitely, and it is more expensive (that is, more costly than O(N)) to change s or 7. 

We have not described the details of how this iff -iff connection problem relates to the two issues presented 
at the beginning of this section (i.e., using the FFT and modification of s for the weighted functions ^^(x)). 
The problem of using the FFT we will postpone until Part II, which describes computational issues. In 
Section |4~2~3| we will describe a method for modification of the decay parameter s, for which the connection 
process described in this section is an integral part. 

4.2.2 The iff-tp Connection Problem 

We now consider the following problem: let / 6 L 2 ([— tt, tt], <D). We assume 7 > and consider two 
expansions: 



16 



The modal coefficients are defined in the following way: 



Jk 



ft = (/^i 7) 

We assume that the modal coefficients for the uppercase (unweighted function) expansion are known and 
that we wish to determine the lowercase modes . From the definitions of the modal coefficients, it is clear 
that we can rewrite the lowercase modes as 



ftp 

Jk 



= U 



,(7) 



(7) 



That is, the modal coefficients for the lowercase basis are identical to modal coefficients of a different function 
for the uppercase basis. To see how this helps us, we make a small digression; recall (pi and define 



9(0) ■■= f 



(7) 



fx 



i(l + e-' w ) 

Suppose that 7 = G e N and that we can somehow find the modal coefficients 



(18) 



-*,(o) 
9 k 



Then we can use the \I r -\f r connection problem outlined in Section [4. 2. 1| to accurately and efficiently determine 
the modal coefficients gf for 7 = G due to the sparse connection. To see how we can find the modal 

Then (18 1 implies that 

G 

(19) 



coefficients , assume that we have the modal coefficients f? 



G 

£ 

m=0 



G 

rn 



.*,(o) _ ;*,(o) 

■Vfe+m — Jk 



If we assume a finite expansion so that g^ = for |fc| > 2iV+ 1, then we can solve ( 19 1 via back-substitution 



Note that determining each coefficient costs 0(G) operations, independent of N; this is a similar operation 
count to the W-W connection cost. 



Finally, we must obtain /*' ( ' ' ) from the given input /*'^- 



However, this is another W-ll/ connection 



(albeit in reverse). Therefore, the three steps to take us from /*'*' G ^ modes to ff'^ modes are 
1. Compute /*' ( ' ^ from J*'^, which is a (backward) connection 



2. Compute g*'^ from /*'^ using (j 1 9 ) 



3. Compute f*>™ = fi>™ 



from g^^°\ a (forward) ^'- v l' connection. 



This is illustrated in Figure [5] For an expansion with N modes, all three steps have 0(NG) cost asymptot- 
ically. The backward connection problem (determining /*>( G ) from /^>( G )) is also computable in 0(NG) 
operations, and is accomplished by reversing the above operations. 



Note that if 7 ^ Nq then all of these steps break down: the connection is not sparse, and ( 19 ) is 



not valid since 7 is not an integer in ( 18 ) 



Performing these modal connections on the real line for expansions in ^ s ^{x) and cj)^(x) is equivalent, 
except one must assign 7 := s — 1 and then proceed as outlined above. 

This particular connection problem is not necessarily useful explicitly since in many of our applications, 
we will have direct access to /*>(°), but each of the pieces necessary for this computation are used extensively 
both in modification of the decay parameter s and application of the FFT. 
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Figure 5: Flowchart representation of a ty-ip connection. The operator -= is the subtraction-assignment operator 
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7 V>,(3) = -*,(3) 
Jk — ilk 





^ - ^ connec- 




tion Figure [4] 



9k — gk+i 



* (0) modes of 
x / 





iff - iff connection 








n k — Jk 


Figure [4] 


► 







= f x \/iu( 3 ) 



h = f x Vw&i 



Figure 6: Flowchart of operations for modification of s. The operator -= is the subtraction-assignment operator. 



4.2.3 Modification Of s: The ip-ip Connection 

We have now developed the necessary tools for the modification of s, i.e., the ip-ip connection problem. We 
assume that G, F £ N and that we know connection coefficients of some function / e L 2 for an expansion 
in ip( F \ and wish to obtain the coefficients for a tfj^ expansion. The whole procedure can be accomplished 
in three steps: 



1. Obtain expansion coefficients for / x 
(iff-ip connection) 



,(G) 



in the * (0) 



in the * (0) 



2. Obtain expansion coefficients for / x 
(Fourier connection) 

3. Obtain the sought expansion coefficients of / in the 
(iff-ip connection) 



Step 2 is easily performed using a version of ( 19 1 by noting the relation between / x 
-n -l 



(G) 



and 



/x 



with knowledge of the canonical Fourier expansion coefficients (iff^(9)). This is shown in 



Figure [6] for the special case F = 3, G = 5. 
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Note that this particular connection problem is very amenable to an FFT+collocation approach whereas 
the algorithm we have laid out is a 'Galcrkin' approach. The problem with the collocation approach is that it 
requires O (NlogN) operations with two FFT's, whereas the above algorithm requires only O (NG) steps. 

As with the *f>-ip connection of the previous section, if either the starting parameter F or the target 
parameter G are not integers, then this procedure cannot be used: the core of the fast algorithm is the 
ability to obtain the canonical Fourier modes, which cannot be done efficiently if the decay parameters are 
not integers. 



4.3 Quadrature 

We now turn to quadrature rules that will compute integrals over the real line. We adopt the following 
notation: the pair jri a,/3 \ denotes the TV-point Gauss-quadrature for the Jacobi polynomial of 

I J 71=1 

class (a, /?), i.e., 

f f(r)w^dr = J2 f (r { n' 0) ) V/ e B 2N . U 

J - 1 71=1 



where B2N-1 is the space of polynomials of degree 27V — 1 or less. We suppress the dependence of r„ 



(«,/?) 



and oj^n'^ on N. We also denote 



the fixed node r 



(a,/3);GR _ 



GR ; ,(a,/3);GR\ 



N 



J n=l 



as the TV-point Gauss-Radau quadrature with 



N 



= 1. We assume for clarity of presentation that the nodes are ordered by n, e.g. 



With the goal that we wish to develop quadrature rules for the infinite line, we will take pains to develop 
quadrature rules in 0-space that do not have nodes at 9 = ±7r, which map to x = ±00. We use the 
Jacobi-Gauss quadrature rules as the building blocks for our generalized Fourier quadrature rules. 

Suppose we wish to construct an TV-point quadrature rule associated with the functions *S>^ (9) . If TV is 
even, then define 



— arccos 



( r (-V2,7-V2)) 



1< n < 



N 



-e 



(7) 

JV+l-n' 



+ 1 < n < N, 



where ri"'^ comes from an ^-point quadrature rule, and 



(20) 



(-l/2, 7 -l/2) 
OJn , 



"JV+l-n' 



1 < n < 



N 



N 



2 + 1 < n < N, 



and Wn*'^ comes from an y -point quadrature rule. 
If N is odd, then define 

^(-l/2,7-l/2)|GRj ^ 



9^ = 



arccos 1 



V Cjv+l-n' 



where r%*'^' GR comes from an ^±i-point quadrature rule, and 



1 <n< 



^<n<N, 



(21) 



(-l/2, 7 -l/2);GR 



2wi" 1/2 ' 7 ~ 1/2);GR , 



1 < n < 



N-l 



_ N+l 



(7) 

N+l-n' 



Tl = 



^<n<N. 
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even 



Jacobi-Gauss Quadrature 




[N odd) 
Jacobi-Gauss-Radau Quadrature 




Figure 7: Construction of Gauss- type quadrature for generalized Fourier functions. The new quadrature rules are symmetric 
combinations of Jacobi-Gauss-type quadrature rules. The constructions shown are accurate node locations for 7 = 5. 



For graphical descriptions of the above formulae, see Figure [7] We have used Jacobi-Gauss rules for N even 

and Jacobi-Gauss-Radau rules for N odd. By construction, when N is odd, 6^+1 = due to the Gauss- 

2 

Radau rule requirement that r^v_[^' — 1. Xhe quadrature rules derived above have no nodes at 9 — zb7r 

2 

(since there are no Jacobi-Gauss, or Jacobi-Gauss-Radau nodes at r = —1) and are symmetric rules for any 
7. Thus they are always exact for any odd function. It is not difficult to show the following result: 



Proposition 3 For N even, the N '-point quadrature rule satisfies 

fix N 

/ e tke w^ 0) d9 = ^ ^\ 

J —7T 1 



x ... , „ \k\<N-l. 

n=l 

When N is odd, the quadrature rule satisfies 

i:„ e^w^ 0) d6 = E^ =1 e™^ (d™) ni 7) , \k\<N. 

The degeneracy in the quadrature rule for N even is exactly of the same nature as the degeneracy in the 
canonical equispaced Fourier quadrature rule for an even number of grid points |18j . If 7 = the rule 

is exactly the same as the equispaced Fourier quadrature rule, symmetric about 9 = 0. The 
quadrature rule {#n°\ f^ "*} can be used to integrate against the weight function Wg 1 ^ when 7 € N 

I J n— 1 

since in this case the weight is itself a trigonometric polynomial. 

In order to determine a quadrature rule to integrate the weighted functions -^^(9), we can augment the 

weights f2n 7 ^ to contain information about the weight function. This can be summed up in the following 
result: 

Corollary 3 The even N -point quadrature rule |#i 7 \wi 7 ^| , where LOn^ := w g 7 ' ' y9„^ ^n 7 '' satisfies 



^^dO = E ^ (W) ^ ^\ 1*1 + \l\ < N - 1 

n=l 

Multiplying fil 7 -* by the inverse of the weight u^ -7 ' ^ is mathematically not a problem since none of 
the 9n^ are equal to ±71", where the weight Wg _7 '°- ) is singular. Note that since the functions $^\x) are 
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s = 7 + 1 



o^°°% o^ a °° 



O 

o 



j 7=0 f I 7=2 f ! 7=4 £ I 7=6 J [ 7=8 ° 

• • • t> o \ » 

b o ^o° U o 00O o°' 



s = 1 < • — — • — 

s = 3 — • 
s = 5 < — o 

s = 7 — 
s = 9 — 

Figure 8: (Top) Plots of the Fourier quadrature nodes on the unit circle generated with equation \2l\ , N = 21. (Bottom) 
The resulting quadrature nodes on the real line. The scale on the real line is \x\ < 15. 



just a mapping of the functions ^^7\9), the quadrature rule , which has nodal 

values over R, can be used to integrate the functions (x) over the real line. Similarly, the rule 

can be used to integrate Galerkin products of the generalized Wiener functions 

cj)^\x) over the real line. 

For various 7/s we graphically depict the location of the quadrature nodes for = 21 in Figure [8] on 
the unit circle z £ T and on the real line. Note that as we increase 7 the quadrature nodes become more 
and more concentrated towards z — 1 (9 — 0). On the real line, this manifests itself as higher concentration 
near x — which, although rcctifiable via an affine mapping, is suboptimal if one wishes to resolve functions 
away from x = 0. Note that the tendency of Jacobi-Gauss nodes to become more equidistant on [—1, 1] as 
(3 (i.e. 7 or s) is increased [18] also suggests that these generalized quadrature rules for large 7 or s will not 
be as good as the the ones for smaller 7 or s since equidistant nodes are bad for finite-interval polynomial 
interpolation. In addition, when 7 = 0, we can use these (equidistant) quadrature nodes to employ the FFT 
for modal-nodal transformations. 



4.4 The Stiffness Matrix 

In many applications to differential equations it is necessary to express the derivative of a basis function as 
a linear combination of basis functions. We devote this section to this endeavor. We define entries of the 
stiffness matrix as 




For the generalized Wiener rational functions, the following result can be proven: 

Theorem 2 Let denote the N x N stiffness matrix for the weighted Wiener rational functions <j>^ . 
satisfies the following properties for any s > ^ •' 
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1. is skew-Hermitian, i.e. Sf! l — —Sf k 

2. is sparse with entries only on the super-, sub-, and main sinister and dexter diagonals: define 

fc v := sgn(fc) (|fc| - 1) = k - sgn(fc), k A := sgn(fc)(|/c| + 1) = k + sgn(fc). 

Then 



ie{±fe v ,±fc,±fe A } 

( s) 

for some purely imaginary constants / . In other words, 



s£, = 0, I g {±fc v ,±fc,±fc A }. 
3. The spectral radius of satisfies 

piS*) <N + 5s. 

The proof of Theorem [2] is quite tedious, so we only sketch the main points. Details are given in Appendix 

m 

Proof Property 1 can easily be deduced by using integration by parts and noting that the functions <j>^ (x) 
decay to zero as \x\ — ► oo. 

Property 2 is a highly nontrivial result that is provable using several properties of Jacobi Polynomials. 
We refer the reader to [24 . Most of the calculations are straightforward once a list of Jacobi Polynomial 
properties has been compiled. However, there are some difficulties whose resolutions rely on a couple of 
fortuitous properties: first, that d x jp (6) = 1 + cos 6*, i.e. that the map we have chosen to take 6 — ► x has a 
Jacobian with a particular form. Second, that 

^[(sin^+W) ( COS 0)' 

is a sparse combination of P^ a '' 3+1 ' ) (cos 9), which is not a trivial result; we show this by using brute-force 
calculation with the compiled list of Jacobi Polynomial properties. 

Property 3 can be derived from the second property. The key ingredient is Gerschgorin's Theorem. Using 
the explicit entries for the constants Tj, { given in Theorem [3] of Appendix [i] we can show that for each k 
satisfying |fc| > 2 the following crude bounds hold: 

|Tfe,fe! < n + 2s, 

|Tfe,_fc| + |rfc,fcv| + |rfe,_fev| + |T fe)fc A| + |r fel _ fc A| < n + 3s + 2, 

where n :— \k\ — 1. Gerschgorin's Theorem can now be used to define a region in the complex plane 
in which all the eigenvalues lie. By the above properties, this region has distance from the origin at most 
In + 5s + 2. Once we consider the necessary relationship between n, k, and N, the result is proven. (It is 
interesting, but not necessary, to note that the eigenvalues all lie on the imaginary axis due to the skew- 
Hcrmitian property of S.) □ 



Remark 2 While the O(N) maximum eigenvalue does depend on s, the proportionality factor is empirically 
around 2, not 5 as given in the theorem. See Table [2j 

The sparsity pattern we have derived for the derivatives of these functions (property 2 of the above 
theorem) is illustrated in Figure [9] Note that the unweighted functions <&( s )(x) also have a similar sparsity 
result; see Lemmakl However, the Fourier functions ^^(6) and tp^{6) do not have sparse stiffness matrices 
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Figure 9: Sparsity plots for stiffness matrices of the weighted Wiener rational functions <pY ' . The sparsity patterns are 
representative of property 2 in Theorem [2] for s = 1 (left) and all s ^ 1 (right). The s = 1 sparsity pattern has been derived 
previously I13| . and the expressions for the r^j in Appendix fBl with s = 1 reduce to the pattern above. 



(unless 7 = 0). In addition, numerical values for the maximum eigenvalues of the stiffness matrix (property 
3) are given in Table [2] The sparsity of the stiffness matrix is important for fast computations of derivatives 
for spectral methods for solving PDEs, and the 0(N) maximum eigenvalue of the stiffness matrix indicates 
that we can take a relatively large timesteps for time-dependent problems. Finally, the skew-symmetry of the 
stiffness matrix easily leads to energy conservation for the Galerkin discretization of hyperbolic conservation 
laws. 



5 The Semi-Infinite Interval 

The generalized Wiener basis functions we have derived can be used for function expansions on the infinite 
line. In order to address expansions on semi-infinite intervals, we can instead use either the even or odd 
Jacobi polynomial basis sets that make up the Fourier functions constructed in Section |2.3| 



The Jacobi functions from Lemma [T] can be mapped and weighted in a procedure identical to the con- 
struction of the Wiener basis. The result is the collection of functions 



8/2 5C-1/2..-3/2) ( ,-. 

X 2 +1 1 r n 1 1+x2 1 , n t « 



(22) 
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s\N 


11 


50 


101 


250 


501 


0.6 


7.31 


43.76 


91.50 


237.60 


483.75 


1.0 


7.99 


44.51 


92.28 


238.39 


484.54 


6.0 


15.96 


53.75 


101.81 


248.14 


494.40 


7T 2 w 9.87 


21.72 


60.67 


109.05 


255.63 


501.99 


15.5 


29.73 


70.45 


119.40 


266.44 


512.99 



( s ) 

Table 2: Maximum eigenvalue of the N X N stiffness matrix S 9 for the Wiener rational functions <p k . The results adhere to 
the asymptotic bound given in property 3 of Theorem [2| 



These functions are a direct mapping and weighting of the Jacobi polynomials. Because of this, they are 
orthonormal and complete in L 2 (]R + ,]R). Mapping techniques for classical functions are not novel and we 
discuss existing methods in Section [6] The classical competitor for spectral expansions on semi-infinite 
intervals is the set of Laguerre functions (weighted Laguerre polynomials). A comparison between the 
Laguerre functions and the functions defined in (22 1 will be made in Part II, and in Section [6] a different 
mapping transformaing Jacobi polynomials to the semi-infinite line will be addressed. 



We make use of the regular square root function 



in (22 1 instead of the phase-shifted version 



,0) 



because there is no need to have complex-valued functions. The phase-shifted square root was a convenient 
choice for the Wiener functions on the infinite line: its compact Fourier series representation ^ enabled fast 
connections (Section 4.2 1 and sparse differentation matrices (Section 4.4 1. By using the real- valued square 



root in (22 1 we sacrifice these two properties. However, the FFT can still be used for the evaluation of modal 



coefficients if s is an integer. 

The caveat in using these functions for expansions on the semi-infinite interval is the fact that they all 
have zero- valued odd derivatives at x = 0. This parallels the same property at 6 = for a cosine series 
on 9 £ [0, it]. Alternative mappings of the Jacobi polynomials to the semi-infinite line do not exhibit this 
restriction, but those mappings also preserve the 0(N 2 ) time-stepping restriction for nodal-based polynomial 
solvers of time-dependent partial differential equations using explicit time-integration on finite intervals. In 



constrast, the functions (22 1 only have an O(N) time-step restriction, similar to the time-step restriction for 
a finite-interval cosine basis expansion. 



The restriction of the Wiener functions to the semi-infinite interval as defined in (22 1 comes both with 



advantages and sacrifices. These functions are purely weighted maps of Jacobi polynomials and are therefore 
easy to implement. Some of the attractive features of the Wiener rational basis functions on infinite intervals 
are lost (e.g. sparse stiffness matrices). However, these functions have properties that are advantageous when 
compared with existing mapping techniques (Section p3J). A numerical comparison between those mapping 
techniques, the functions (22 1, and the Laguerre functions will be made in Part II. 



6 Alternative Methods 

Before concluding this article with a summary of the derived properties of the generalized Wiener basis, we 
first summarize existing results on the topic of mapping Jacobi polynomials from the finite interval to the 
infinite interval. This method is very closely related to our strategy of mapping a generalized Fourier series 
from the canonical finite Fourier interval to the real line. Numerical studies comparing these methods are 
presented in Part II, but it is appropriate to acknowledge these functions here, and to discuss how they relate 
to the Wiener rational function basis. 
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6.1 The Infinite Interval 



The main idea for our generalization of Wiener's original rational basis is using a 'well-behaved' mapping 
to transform functions on a finite interval to those on an infinite interval. This basic idea is classical |15j . 
Indeed one of the more popular mappings that has gained momentum in the literature are the so-called 
'mapped Chebyshev' functions/polynomials. 

In order to further generalize the mapped Chebyshev functions, we will briefly restate their derivation 
from our point of view. We begin with the Jacobi polynomials Pn° (r) on r G [—1,1]. Mapping via 
r = cos 8 to 8 £ [0,7r] yields trigonometric polynomials. We now 'stretch' the domain to € [— 7r,7r] via the 
affinc mapping O = 28 — it. Finally we utilize the usual linear fractional map e 10 = (i.e. rotation of 
the Riemann Sphere) to yield functions on the real line i6R. For all s,t > |, this results in the functions 
PB^'^O), defined as 



(x) = P ? (( 2;s - 3 )/2^2t-3)/2) { 

orthonormal on the real line under the weight 



(«,*) 

' J PB 



-i (2s-3)/2 r 



1 - 



l 



(2t-3)/2 



and the weighted functions 



(*,*) ._ 



are orthonormal under the unweighted inner product. When s = t = 1, the functions PB^' f ' coincide 
with the mapped Chebyshev polynomials T~B n (x) introduced in [3] and subsequently developed in [5] and 
9 , although the original idea of applying spectral expansions over finite intervals to solving problems over 
infinite intervals seems to come from [16] . In any case, the mapped Jacobi functions pb^ 8 '*^ decay like y^p 
for x — > — oo and tAj for x — ► +oo. The advantage of these functions is that the decay can be different as 
\x\ — > oo. Also, others have already explored some convergence theory in function spaces [3] and applications 
to differential equations |31j for the Chebyshev case s = t = 1. In Part II when we present numerical 
examples, we will use the basis set pb^ 5 ''-* with s = t = 1, i.e. the Chebyshev case. 

Note that because all of these mapped types of polynomials and the generalized Wiener basis we have 
presented ultimately stem from Jacobi polynomials and mappings of similar character, all these basis sets are 
related in some fashion. To relate the mapped Jacobi functions to the generalized Wiener rational functions, 
we have 

PB(-)(x)«He(*« 

\x - \JX A + 1 + 



In Table[3]we relate the unweighted functions to the generalized Wiener rational basis, modulo multiplicative 
constants. In this article we make no observations about how mapped Jacobi polynomials compare to the 
Wiener basis set as a practical tool for function expansions. However, such a comparison will be a central 
theme in Part II. 



6.2 The Semi-Infinite Interval 

To perform spectral expansions on semi-infinite intervals, the only classical technique is the Laguerre poly- 
nomial/function method. However, mapping techniques can be used to transform finite-interval methods to 
semi-infinite interval methods. 

As with Section [6. 1[ we explain the choice of mapping from our point of view as a mapping of the Riemann 
Sphere. The Jacobi polynomials are defined on r G [—1, 1]. If we allow complex values of r, then we may 
consider using a linear fractional map to transform the Jacobi polynomial domain to the semi-infinite line. 
The ordered assignments r = {1, 0, —1} to x — {0, 1, oo} specify the transformation uniquely as 
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Previous function 


Name / classification 


Interval 


Reference 


Relation 




Chcyshev rational functions (1st) 


R 


uni, H, El 


pb! 1 ' 1 ) 

n 


& / U JDti 
it/ 1 1 


Chebyshev rational functions (2nd) 


R 


ITCH, El, M 

11 ' (-H ' MH 


p B (2,2) 

n 




Christov functions (even) 


R, [0, oo) 


[TBI. HTJ1 

1 — -—i * 1 — —4 


Im{^ 0) } 


*~ r 1 i / 1 i 


Christov functions (odd) 


R, [0, oo) 


1151. fTTjl 




CH„ 


Higgins functions (even) 


R, [0, oo) 


□n 

1 — --H 


Rcj^ 1 ' '} 


SH n 


Higgins functions (odd) 


R, [0,oo) 


pm 


Im{^ 0) ) 
1 J 


Pk 


(Complex) Higgins functions 


R 


m, m 




Cjfc 


(Complex) Wiener rational functions 


R 






TL n 


Half-infinte Chebyshev rational functions 


[0,oo) 


m 


p L (l/2) 



Table 3: Relationship between orthogonal functions in previous work and the current bases presented. 



U r=±§. (23) 



If necessary, one can also specify the relationship to 9 and the cosine series on [0, %]. For details, see 0. Our 
definition of the transformation differs only in orientiation from that presented in [7] . We have chosen this 
orientation so that the Jacobi parameter (3 is assigned to the location x = oo in order to mimic to the same 
assignment for the Wiener functions. 



In the literature the maps of the Chebyshev polynomials under the transformation ( 23 1 are labeled 
TL n (a;). Adopting similar notation, we define 



PL W (x) = i*-W*-*) 



G [0,( 



which are L 2 -complete and orthonormal under the weight function 

(2s) 



It is then possible to define the weighted functions 

pl - )(x) = (rr^) SpL " )(x) (24) 



2 ' S ^- lA2s - 2) f^) 5 (25) 



1 + x ) n \l + x 

which are L 2 -complete and orthonormal under the weighted L 2 inner product 

for any s > 4. The pl^ s ' are defined for x <E [0,oo] and decay like x s as x — ► oo. A significant difference 
between the Wiener-type functions (both on the infinite and scmi-infinte intervals) and the set defined in 
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(24 1 is the fact that these functions are not orthogonal in the unweighted I? inner product, but instead in 
the norm defined by the above inner product. This choice was made (as opposed to defining functions in 
the unweighted inner product) to ensure that integer values of s resulted in a Jacobi polynomial family that 
was amenable to usage of the FFT. 

The main observation we make regarding this basis is that these functions are the result of a linear 
fractional map directly from the Jacobi domain; therefore, they will inherit the 0{N 2 ) time-step restric- 
tion of nodal explicit time-integration methods for time-dependent partial differential equations. The same 
observation can be made about the functions defined in [7]. 

7 Conclusion 

We have presented a collection of generalized Fourier series which, when mapped and weighted appropriately, 
generates a basis set on the infinite interval with a tunable rate of decay. For each rate of decay s satisfying 
s > \ the resulting basis set 4>k '■ 

is orthonormal and complete in L 2 (R, <D) 
is characterized by x~ s decay for \x\ — ► oo 

• can be generated via Jacobi polynomial recurrence relations 

• has sparse connection properties that can be efficiently exploited via combinations of sparse Fourier 
and Jacobi connections 

• has an N x N Galerkin stiffness/differentiation matrix that has at most 6 A nonzero entries with O(N) 
spectral radius 

• is characterized by a 'Gauss-like' quadrature rule. 



When s € N, the basis set is a rational function; we will show in Part II that in this case we can use the FFT 
for modal-nodal transformations. The case s = 1 corresponds to a mapping and weighting of the canonical 
Fourier series, as discovered by others previously. Due to the original presentation of the s = 1 basis by 
Wiener [3U], we call the functions tfi^ the generalized Wiener rational basis functions. 

These basis functions have a similar flavor to directly mapped and weighted Jacobi polynomials (called 
pb£ s,t) here). In Part II we will compare these basis sets and discuss advantages and disadvantages of each. 
In addition, we will also employ the Sine and Hcrmite functions in test cases in an attempt to investigate 
a relatively broad class of spectral approximation methods. In contrast to |26j which reviews much of 
the theory present for expansions on the infinite interval, we concentrate on numerical issues, including 
application of the FFT. We will extend our investigation to the semi-infinite interval to compare the Laguerre 
polynomials/functions, the mapped Jacobi functions (denoted pl„ here), and the restriction of the Wiener 
functions to the semi-infinite interval as given in Section [5] 

We do not wish to claim that, on the infinite or semi-infinite intervals, genuinely global spectral expansions 
are truly superior to alternative numerical approximations; rather we wish to identify the generalized Wiener 
basis set as a novel competitor to existing global spectral expansions. Part II will follow up to show that the 
Wiener basis set is very competitive with existing expansions. 
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A Recurrence Coefficients 



In this appendix wc compile various recurrence relations for the Jacobi/Szego-Fouricr/Wiener rational func- 
tions. We state the recurrences in terms of the Szego-Fourier functions ^ k , but note that they all apply 
equally well to the unweighted Wiener rational functions as well. Note that we only list recurrences for 
k > 0; for k < 0, we may use the conjugation relation jlo| ) to obtain V?^^ a t almost no additional cost. We 
first require a tour of some Jacobi polynomials recurrences: 



,(a,/3) fi(a,/3) 
u n+l r n+\ 



* n 

(1 - r)Pt'V 

± n 



dr 



p(a,f3) _ J h (a,f3) fi(a,/3) 
V n n— 1 ' 



(a,/3) p( Q -l,/3-l) 



i=0 
2 



(a,/3) 5(a+l,/3+l) 



n. — i n—i 



„(«,/3) p(a-l,/3) _ „(<*,/?) p (a-l,/3) 



,,{fta)pM-l) , „(/3,«) p(«,0-l) 
H>n,Q r n T Mn,l r n+l ' 

,.(«,/?) p(a+l./3) _ p(a+l,/3) 
„(/9.«) p(a./3+l) , ,.(/9,a) p(a,/8+l) 

(a,/3) p(a+l,/3+l) 

m J n — 1 1 



(26) 
(27) 

(28) 
(29) 
(30) 

(31) 
(32) 

(33) 



where ^n o/v ^n°o/-l> an< ^ 1^'^ m (29 1- (33) are constants for which we take explicit formulae from |25j : 



,K/3) _ 



2(n + a)(n + a + /3) 
(2n + a + /3)(2n + a + /3+ 1)' 

f 2(n + l)(n + /?+!) 

(2n + a + + l)(2n + a + /? + 2) ' 

f 2(n + a + l)(n + a + /3+l) 
(2n + a + /3 + l)(2n + a + /? + 2) ' 

' 2n(n + (3) 

(2n + a + (3)(2n + a + f3+ 1)' 



(34) 

(35) 

(36) 

(37) 
(38) 
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The three-term recurrence coefficients in ( 26 1 are given by [14] : 



P~ a n = 

P 2 -* 2 „ > 

(2n+a+/3)(2n+a+/3+2) ' " ^ U- 

2°+ g+1 r(q+l)r(/3+l) 
r(a+/3+2) ' 

4(q+l)(j3+l) 
(a+/3+2) 2 (a+/3+3) ' 

4n(n+a)(n+l3)(n+a+f3) 

(2n+a+/3-l)(2n+ct+/3) 2 (2n+Q+/3+l) ' 



(39) 



71 = 0, 

n = 1, 
n > 1. 



(40) 



The demotion recurrence coefficients in ( 27 1 can be obtained by determining the analogous relations for the 



monic orthogonal polynomials (pQ, [27]) and then employing the appropriate normalizations: 



-(«,/3) 
-n.O 



q/3 



(a+/3)(a+/3+l)' 



2 

(Q+/3+2) V (a+[3+3) ' 



' (rt+a)(«+ / 3)(n+q+ ) 3-l)(n+a+/3) 
(2n+a+/3)V (2n+a+/3-l)(2n+a+/3+l) ' 



2(<*-ff) 



(a+/3+2) v / ^ : ha' 

2(a-/3) v / (n+l)( ra +a-|-/3) 
(2n+a+/3)(2n+Q+/3+2) ' 

2 



n = 0, 
n > 0. 



2(a+l)Q3+l) 



a+/3+2V (a+/3+l)(a+/3+3) ' 



2n+a+/3+2 



Finally the promotion relation (28 1 coefficients can also be determined: 



'ln,0 



e 



(a+1,/3+1) 

Tl,0 1 



n = 0, 



n = 1, 



n > 1. 



(rt+l)(n+2)(7t+q+l)(n+/3+l) 
(2n+a+0+l)(2n+a+f3+3) ' 



n = 0, 
n > 0. 



(41) 



(42) 



(43) 



'ln,-l 



e 



(a+1,/3+1) 
n-1,1 ' 



(a,/3) _ _ (a+1,/3+1) 
'/n,-2 — fc n-2,2 



Of course, ( 27 l-(|28|are consequences of combining ( 29 1-( 32 1 . Using the orthogonal polynomial three-term re- 
currence relation (|26| we can show the following recurrence relation for the Szego-Fourier functions ^^' S \9): 



(7) 
n+1 



U^cosO-V^ 



^ n 



(7) 



VV n y„_i "-n T -( n -l)- 

In the following expressions, we make use of the following definitions: for a given 7 > — A, 

a:=-\, /3:=7-|. 



(44) 
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The recurrence constants are then given by 



rr(7) 
U ±n 



± 



(c + 1,/3 + 1) 



V. 



(7) 



±n 



± 7 



(«,£) 







(7) 



Using the promotion and demotion three-term recurrences ( 27p8 ) we also have the following recurrence 
relation: 



(7) 

n+l 



where the recurrence constants are given by 



(7) 



(45) 



r~r(7) 
U ±n 



1 =r I 

„(<*■« T (a+l,jB+l) 
''ti,0 fc Ti-l,2 



_(a + l,/3 + l) (a,f>) 
-ti-1.1 I '/n , - 1 

_(a, + l,/5 + l) 31 _(«,« 



ir 



(7) 

±n 



_(c + l,/3 + l) 



-n-1,2 



± 



'n,-2 



Finally, putting these last two recurrences together yields 



(7),t,(7) _ 



n+l 



A 



D 



(7) 



(7) 



.4 



(7) 



— B 



(7) 



(7) 



with the following values for the recurrence coefficients: 



D 



(7) _ 



.4 



(7) _ 

±n — 



2s 



(o./9) 



n = 0, 
n > 0. 
n = 0, 
n > 0. 



(7),t,(7) 



n-1 



^r(7)^r,(7) 



B 



(7) 

±n 



(7) 

±71 



27^2 

v / 7+T' 



2y'(™+ 1 )(«+7-i) 



( 7 4 7 ) + [ 2 V^H^)-l] Ag) 



r 0, 



(27+1)7 

4 T) (7+i)V 7+2 ' 

7 £ n ,o 

. V(»+7-2)(n+7-l)A™ 1 



n = 0, 
n > 0. 
n = 0, 
n = 1, 

n > 1. 
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B The Stiffness Matrix 



We assume the decay parameter s > h is given and we derive a and (3 from the value 7 := s — 1 as in 
Appendix [X] Also, we define increments and decrements of the integer index k £ 7L: 

(X . 21 P ■ s 2' 

fc v = sgn(fc) (|fc| - 1) , fc A = sgn(fc) (|*| + 1) , 

n := \k\ - 1. 

(s) 1 

We begin by noting the sparse representation of the product of <pl and 7-3 



(x-i) ■ 



Lemma 3 We have the representation: 



(x-i) 



for some constants Xu i ■ 
Proof We first note that 



-- [sin0(x) + i(l + cos9(x))} 



after making the transformation to 6{x). Then making the identification <&j^ = ^j! 8 , we may use recur- 



rence relations ( 44 1-( 45 1 to obtain the result 



□ 



A second more potent result is the sparsity result for the unweighted Wiener rational functions <&^\x): 
Lemma 4 We have 



d®[ s) (x) 
dx 



l^le{±k v ,±k,±k^} a k,l ' 



with 



J k,±k^ 



?s s n ( fc )2(irq 



+s) y (2n+s-l)(2n+s+l) 



y/(n + s - l)(n + s) ± x /ra(ra + 1) 



x (s) 



i sgn(fc) 



v/(n+ l)(n + .s), 
V(n + l)(n + a) (2ra+ *) ( ( 2 ~ S j s+ , ; ■ 



+k, 
-k. 



1 n+l 



Proof This result can be proven by brute- force calculation of the derivatives using ( 29 )-( 33 1 and the 



recurrence formula ( 26 1 . Two critical steps are necessary: a highly nontrivial collapsing of a special arithmetic 



combination involving various constants in several Jacobi polynomial relations, and the very special form of 
the 8 — > x Jacobian for the mapping. Thus, the particular form of the mapping is critical in proving this 
result. □ 

Putting the two lemmas together, we have the desired sparsity result for the (f>^\x) stiffness matrix: 

Theorem 3 The following equality holds for any s > | : 



dx 



E 



T k,l Wl ) 



ie{±fe v ,±fe,±fc A } 
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(s) 

where the constants ri / are equal to 



Tk,±k v — 



f 1 s(s-2) 

4\/ (2n+s-l)(2n+s+l) 



jsgn(fe) (V( n + s - + s) ± V n ( n + !)) + 
^ (v^ + l)(n + s - 1) ± ^ + i)) } , 



Tfc.fe 



(2 

? sgn(fc) v /(n + l)(n + i) - 2(2n "^_ 

is{s—l) 
2(2n+s)(2n+s+2) ' 



2(2n+s)(2ra+s+2) 2 ' 



M s(*-2) y 

4 V x (2n+s + l)(2n+s+3) 



{~^+2 [ V(« + 2)(n + s) ± v^n + l)(n + s + 1) 
sgn(fc) y/{n + l)(n + 2) ± ^/(n + s)(n + s + 1) | . 



Clearly for \k\ — I we have Tk,±k v — Tfc,o so that, taking into account the different normalization constant in 
the definition of <fJj* for \k\ = 1 we have: 



Tfe ' o = ivlrl {sgn(fcV ^ 1} ' 



ifci = i. 



And for k = 0: 



da; 



2\ s- -. 



Proof The hard work was completed in Lemma [4] The result of this theorem follows from a simple 
application of the product rule of differentiation to 



2 s/2 



dx \_(x — i) s 



with the appropriate use of Lemmas [3] and [4] 



□ 



With explicit entries for the stiffness matrix derived, we are now ready to give the proof the third property 
of Theorem [2] the spectral radius of the stiffness matrix. 



33 



Proof We can crudely bound the entries of the stiffness matrix for n := — 1 > 1 





< 


f +S, 




< 


S 

2 ' 


Tfe,±fev 


< 


f + s, 


T~k,k\ 


< 


n + 2s, 


Tk,-k\ 


< 


1, 


Tk.kA 


< 


f + S + 


Tk,-k A \ 


< 


S 

4' 



Thus we have: 

|T_ fc | + |r fe v| + |T_ fe v| + |r fe A| + |r_ fc .| < 1 + f + s+f + f + s+i + f 

= n + 3s + 2. 

An application of Gerschgorin's Theorem then proves the result. □ 
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